﻿using System;
using System.Text;

namespace BSU.Math.Algebra
{
    public class Matrix
    {
        #region Private fieds

        private double[][] _a;

        #endregion

        #region Properties

        public int M { get; protected set; }
        public int N { get; protected set; }

        #endregion

        #region Constructors and initializers

        protected Matrix()
        {
        }

        public Matrix(int m, int n)
        {
            M = m;
            N = n;

            _a = new double[M][];
            for (var i = 0; i < M; i++)
            {
                _a[i] = new double[N];
            }
        }

        public Matrix(Matrix a)
            : this(a.M, a.N)
        {
            for (var i = 0; i < a.M; i++)
            {
                for (var j = 0; j < a.N; j++)
                {
                    this[i, j] = a[i, j];
                }
            }
        }

        #endregion

        #region Operators

        public double this[int i, int j]
        {
            get
            {
                return _a[i][j];
            }
            set
            {
                _a[i][j] = value;
            }
        }

        public static Matrix operator +(Matrix a, Matrix b)
        {
            if (a.M != b.M || a.N != b.N)
            {
                throw new ArgumentException("Cannot perform operation with matrices with given size.");
            }

            var c = new Matrix(a.M, a.N);

            for (var i = 0; i < a.M; i++)
            {
                for (var j = 0; j < a.N; j++)
                {
                    c[i, j] = a[i, j] + b[i, j];
                }
            }

            return c;
        }

        public static Matrix operator -(Matrix a, Matrix b)
        {
            if (a.M != b.M || a.N != b.N)
            {
                throw new ArgumentException("Cannot perform operation with matrices with given size.");
            }

            var c = new Matrix(a.M, a.N);

            for (var i = 0; i < a.M; i++)
            {
                for (var j = 0; j < a.N; j++)
                {
                    c[i, j] = a[i, j] - b[i, j];
                }
            }

            return c;
        }

        public static Matrix operator *(Matrix a, double k)
        {
            var c = new Matrix(a.M, a.N);

            for (var i = 0; i < a.M; i++)
            {
                for (var j = 0; j < a.N; j++)
                {
                    c[i, j] = a[i, j] * k;
                }
            }

            return c;
        }

        public static Matrix operator *(double k, Matrix a)
        {
            var c = new Matrix(a.M, a.N);

            for (var i = 0; i < a.M; i++)
            {
                for (var j = 0; j < a.N; j++)
                {
                    c[i, j] = a[i, j] * k;
                }
            }

            return c;
        }

        public static Matrix operator *(Matrix a, Matrix b)
        {
            if (a.N != b.M)
            {
                throw new ArgumentException("Cannot perform operation with matrices with given size.");
            }

            var c = new Matrix(a.M, b.N);

            for (var i = 0; i < a.M; i++)
            {
                for (var j = 0; j < b.N; j++)
                {
                    for (var k = 0; k < a.N; k++)
                    {
                        c[i, j] += a[i, k] * b[k, j];
                    }
                }
            }

            return c;
        }

        public static Vector operator *(Matrix a, Vector b)
        {
            if (a.N != b.M)
            {
                throw new ArgumentException("Cannot perform operation with matrix and vector with given size.");
            }

            var c = new Vector(a.M);

            for (var i = 0; i < a.M; i++)
            {
                for (var j = 0; j < b.M; j++)
                {
                    c[i] += a[i, j] * b[j];
                }
            }

            return c;
        }

        #endregion

        #region Public methods

        public override string ToString()
        {
            var s = new StringBuilder();

            for (var i = 0; i < M; i++)
            {
                for (var j = 0; j < N; j++)
                {
                    s.Append(this[i, j]);
                    if (j != N - 1)
                    {
                        s.Append("\t");
                    }
                }
                if (i != M - 1)
                {
                    s.Append("\n");
                }
            }

            return s.ToString();
        }

        public void SwapRows(int i, int j)
        {
            if (i < 0 || i > N - 1 || j < 0 || j > N - 1)
            {
                throw new ArgumentException("Row index is not valid.");
            }

            for (int k = 0; k < M; k++)
            {
                var tmp = this[k, i];
                this[k, i] = this[k, j];
                this[k, j] = tmp;
            }
        }

        public void SwapColumns(int i, int j)
        {
            if (i < 0 || i > M - 1 || j < 0 || j > M - 1)
            {
                throw new ArgumentException("Column index is not valid.");
            }

            for (var k = 0; k < M; k++)
            {
                var tmp = this[i, k];
                this[i, k] = this[j, k];
                this[j, k] = tmp;
            }
        }

        public Matrix Transpose()
        {
            var c = new Matrix(N, M);

            for (var i = 0; i < M; i++)
            {
                for (var j = 0; j < N; j++)
                {
                    c[j, i] = this[i, j];
                }
            }

            return c;
        }

        public Matrix Inverse()
        {
            var n = M;

            var x = new Matrix(n, n);
            var p = new Matrix(n, n);

            Matrix c = LUPFactorization(this, out p);

            for (var k = n - 1; k >= 0; k--)
            {
                x[k, k] = 1;
                for (var j = n - 1; j > k; j--)
                {
                    x[k, k] -= c[k, j] * x[j, k];
                }
                x[k, k] /= c[k, k];
                for (var i = k - 1; i >= 0; i--)
                {
                    for (var j = n - 1; j > i; j--)
                    {
                        x[i, k] -= c[i, j] * x[j, k];
                        x[k, i] -= c[j, i] * x[k, j];
                    }
                    x[i, k] /= c[i, i];
                }
            }
            x = p * x;

            return x;
        }

        #endregion

        #region Private methods

        private Matrix LUPFactorization(Matrix a, out Matrix p)
        {
            int n = a.M;

            var c = new Matrix(a);
            p = new IdentityMatrix(n);

            for (var i = 0; i < n; i++)
            {
                double pivotValue = 0;
                int pivot = -1;
                for (var row = i; row < n; row++)
                {
                    if (System.Math.Abs(c[row, i]) > pivotValue)
                    {
                        pivotValue = System.Math.Abs(c[row, i]);
                        pivot = row;
                    }
                }
                if (pivotValue == 0)
                {
                    throw new ArgumentException("Matrix is singular.");
                }

                p.SwapRows(pivot, i);
                c.SwapRows(pivot, i);
                for (var j = i + 1; j < n; j++)
                {
                    c[j, i] /= c[i, i];
                    for (var k = i + 1; k < n; k++)
                    {
                        c[j, k] -= c[j, i] * c[i, k];
                    }
                }
            }

            return c;
        }

        #endregion
    }
}